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Custom-Made  Pyramids 


Shmuel  Peleg 
Oma  Federbusch 
Robert  Hummel 


Abstract 

Pyramids  are  data  structures  used  to  store  and  process  images  at 
multiple  levels  of  resolution.  The  bottom  level  of  a  pyramid  is 
used  to  represent  data  at  a  fine  level  of  resolution,  while  higher 
levels  of  the  pyramid  are  used  for  data  stored  at  coarser  levels 
of  resolution.  For  example,  in  the  Gaussian  pyramid  data  struc¬ 
ture,  each  successive  level  is  obtained  by  local  averaging  and 
subsampling  of  the  immediately  lower  level  in  the  pyramid.  In 
nearly  all  pyramid  implementations  to  date,  the  size  reduction  in 
each  dimension  between  levels  of  the  pyramid  is  a  constant  fac¬ 
tor  of  two. 

This  paper  describes  a  scheme  that  permits  construction  of  py¬ 
ramids  with  arbitrary  size  reductions  between  levels.  The  reduc¬ 
tion  factors  can  be  different  in  each  dimension,  and  differ 
between  levels,  to  adapt  to  a  given  application.  The  user  can 
thus  specify  a  sequence  of  decreasing  rectangular  image  sizes, 
and  construct  pyramids  conforming  to  those  sizes.  Further,  the 
reduction  factors  can  be  made  adaptive  to  region  properties,  ena¬ 
bling  smooth  regions  to  be  reduced  more  than  “busy "regions. 


1.  Pyramids 

Pyramid  data  structures  have  proven  useful  in  the  development  of  image  pro¬ 
cessing  methods  dealing  with  image  features  at  varying  scales  of  resolution.  A  sur¬ 
vey  on  pyramid  structures  may  be  found  in  [1].  Typical  image  pyramids  are  formed 
from  rectangular  grid  arrays  whose  sides  have  lengths  that  are  powers  of  two  in 
extent.  Thus  the  base  level  might  be  512  by  512,  the  next  level  up  would  then  be 
256  by  256,  and  each  successive  level  reduces  each  dimension  by  a  factor  of  two.  In 
the  Gaussian  pyramid  data  structure,  the  reduction  from  one  level  to  the  next  is 
accomplished  by  blurring  the  lower  level  (by  means  of  convolution  with  a  non¬ 
negative  kernel)  followed  by  a  subsampling  of  every  other  pixel  on  every  other  row. 
In  the  Burt  form  of  the  Gaussian  pyramid  [2],  each  level  has  size  of  the  form  2"  +  l 
by  2"  + 1 ,  so  that  in  the  subsampling  operation  left  and  right  edge  pixels  are  included 
on  every  second  row,  and  that  both  the  top  and  bottom  rows  are  sampled.  Nonethe¬ 
less,  the  sampling  rate  is  still  two,  and  we  say  that  the  size  ratio  between  levels  is 
two  in  each  dimension. 

By  subsampling  every  other  pixel  on  every  row  and  offsetting  the  samples  by 
one  pixel  on  successive  rows,  it  is  possible  to  achieve  an  effective  sampling  ratio  of 
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VT  in  each  dimension  [3].  In  this  pyramid  scheme,  every  other  level  is  in  essence 
located  on  a  45°  grid,  and  care  must  be  exercised  when  designing  kernels  and  other 
local  operators  on  these  levels.  Nonetheless,  the  reduced  reduction  factor  between 
levels  results  in  a  better  chance  of  capturing  a  salient  feature  at  an  appropriate  scale 
of  resolution,  at  the  cost  of  doubling  the  number  of  levels. 

In  this  paper,  we  describe  a  method  for  constructing  pyramids  of  arbitrary  size 
at  each  level,  so  that  the  resolution  can  be  reduced  as  needed,  and  not  only  by  fixed 
ratios.  Further,  the  reduction  does  not  have  to  be  uniform  over  the  whole  image, 
and  can  vary  based  on  local  image  properties. 

The  basic  step  in  the  construction  of  the  proposed  pyramid  scheme  is  a  spatial 
resampling  technique  used  in  graphics  [4,5],  related  to  anti-aliasing.  We  will  first 
describe  the  resampling  idea  in  one  dimension.  The  transformation  can  be  applied 
to  rectangular  grids  by  first  sampling  the  rows,  and  then  sampling  the  columns  of 
the  result.  We  will  then  present  a  formulation  of  the  sampling  idea  that  permits  the 
construction  of  pyramids  using  arbitrary  placements  of  pixels  (such  as  hexagonal 
grids).  The  central  property  of  the  sampling  method  is  that  each  pixel  contributes 
fully  to  the  output  samples,  thereby  minimizing  sampling  artifacts. 

As  an  introduction  to  the  resampling  methods  to  be  defined  in  later  sections, 
consider  the  original  sample  points  as  “producers,”  and  the  new  sample  points  that 
form  the  resampled  data  as  the  “consumers.”  The  producers  and  consumers  have 
associated  positions,  corresponding  to  the  sample  point  locations  in  the  domain. 
The  amount  “produced”  by  each  producer  is  pre-specified,  as  is  the  consumption 
level  of  each  consumer.  The  consumers  obtain  their  products  from  a  linear  sum  of 
the  producers,  taking  contributions  of  fixed  amounts,  such  that  the  total  of  all  con¬ 
tributions  by  a  given  producer  is  the  total  amount  produced  at  that  site.  The 
“error”  of  a  sampling  can  be  defined  as  the  sum  over  all  producer-consumer  pairs 
of  the  distance  between  the  corresponding  points  weighted  by  the  contribution  of  the 
producer  to  the  consumer.  The  error  is  minimized  when  the  consumers  use  the 
closest  suppliers  to  compute  their  values.  The  methods  described  below  have  cer¬ 
tain  optimality  properties  with  respect  to  these  measures,  but  we  will  not  pursue  the 
variational  derivation  any  further. 


2.  One-Dimensional  Resampling 

We  first  consider  the  problem  of  resampling  a  vector  (vo,  •  •  •  ,vjv_i)  of  N  pix¬ 
els  to  a  vector  (wo,  *  •  •  ,wm- l)  of  M  pixels.  We  treat  both  cases  M  <  N  and 
M  2=  N. 

2.1.  Uniform  Resampling 

We  use  the  formula 
Ki+n/pl-i 
=  2,  njvj 

j“  l</pl 


where 


.  ITO..I.U  w«u  HU  >3lWCT:r>JW?  FJ'>7y??T^ 
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and 


ry 


min  {  1,  p,  (  p  (y  +  l)  -  i)  }  if  Li/pJ  s  <  i/p 
min  {  1,  p,  (  (i  +  l)  -  P7)  }  if  i/p  s  ;  s  f(i+l)/pl 


To  interpret  this  formula,  first  assume  that  M  <  W,  s o  that  p  is  small.  Then 
nearly  all  ry  for  j  in  the  appropriate  range  will  have  the  value  p,  with  the  exception 
of  the  two  extreme  ry’ s. 

More  generally,  we  can  regard  the  output  range  [0,M]  split  up  into  M  subinter¬ 
vals  [i,  i  +  l],  with  i  =  0,1,  •  •  •  ,  M~\.  The  same  range  is  also  split  into  N  subin¬ 
tervals  [p7,  p  (;  +  l)],  with  j  *0,1,  •  •  •  ,  N  —  1.  The  coefficient  ry  is  the  total 
length  of  the  portion  of  the  interval  [p-y,  p-(./  +  l)]  that  intersects  with  the  interval 
U,  i  +  l]. 

For  the  situation  with  N  *  4  and  M  =  3,  to  resample  (vo,vi,V2,V3)  with 
(wo,w,i,m'2)  we  have  the  formulas 
3  I  1 

wo  =  jvo  +  jvi 

1  j.  1 

»1  =  JV1  +  Jv2 
1  ,  3 

W2  "  JV2  +  — V3  . 


as  shown  in  Figure  la.  Figure  lb  shows  how  the  coefficients  can  be  computed  from 
the  overlapping  intervals. 

Properties  of  this  sampling  method  include: 


N—  1 

•  The  total  contribution  of  any  given  vj  to  all  w’s  is  p;  i.e.,  2)  ry  =  p. 

1=0 
M- 1 

•  The  sum  of  all  contributions  to  any  given  w,-  is  one;  i.e.,  2)  rU  =  1- 

;=0 


»o  h>i  W2 


Figure  la.  Uniform  resampling  weights  for  a  4-vector  to  a  3-vector. 
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Figure  lb.  Method  for  computing  the  uniform  resampling  weights  for  a  4-vector  to 
a  3-vector,  using  the  overlap  of  the  intervals  at  different  resolutions. 


2.2.  Weighted  Resampling 

In  the  uniform  sampling  given  in  the  previous  section,  all  input  pixels  make  the 
same  contribution  p  to  output  pixels.  We  now  extend  the  resampling  to  permit  each 
input  pixel  to  have  an  adaptable  “forward  weight,”  so  that  pixel  j,  with  value  Vj, 
contributes  a  total  of  y.j  to  the  w’s.  These  weight  should  satisfy  the  normalization 
condition 

N- 1 

2  V-j  =  M  , 
j- o 

reflecting  the  desire  to  have  each  of  the  M  output  pixels  to  receive  unity  in  contribu¬ 
tions  from  the  v’s.  When  for  all  j,  we  have  the  uniform  sampling  of  the 

last  section. 


We  use  the  same  principle  as  in  the  previous  section  to  develop  a  linear  resam¬ 
pling  formula 


N— 1_ 

w,  =  2  ruvj  • 
)‘  0 


i- 1  J 

but  now  rij  represents  the  length  of  that  portion  of  the  interval  I  2  I1*-  2  M-t  that 

|*=o  *-o 

intersects  the  interval  [i,  i  +  l].  Note  that  the  input  interval  number  j  has  length  p.;, 
and  that  the  input  intervals  subdivide  the  output  range  [0,A/]. 


We  omit  the  precise  formulas  for  the  r</s,  but  depict  the  situation  for  N  =  4, 
M  —  3,  and  |ios|ii  =  1,  p,2  =  =  0.5,  in  Figure  2.  In  this  case  we  obtain  the 

formulas 


w0  =  v0 

H’  1  =  V] 

1  X  1 

W2  =  yv2  +  yv3  . 


We  can  formulate  the  general  weighted  resampling  formulas  by  giving  an  inter¬ 
polation  formula  and  a  sampling  formula.  Specifically,  given  a  vector 
(vq.vi,  •  •  •  ,vs-i)  and  the  weights  (p,o,  •  ■  •  ,fiAr-i),  we  form  an  interpolation 
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Wi  W2 


Figure  2a.  Pyramid  representations  of  the  weighted  sampling  with  weights  as  in  the 
text. 


M  *  1 


*1  “  1 


(12  =  0.5  =  0-5 


Figure  2b.  Weighted  resampling  of  a  4-vector  to  a  3-vector. 


function 


where 


/(*)  =  2  v<4>«(*)  , 

i“0 


i  —  1  i 

=  •  1  if  2  P-*  s  *  <  2  M>* 

*=0  *  =  0 

‘  0  otherwise  . 

The  samples  (wo,  •  •  •  ,w^-i)  are  then  obtained  from  the  sampling  formula 
w«  =  /  fixWMdx  , 

where 


...  /l  if  i  s  x  <  i  +  l 
[p  otherwise  . 


This  formulation  suggests  a  generalization,  to  which  we  return  in  the  next  section. 
For  the  moment,  we  note  that  the  ?</ s  satisfy 


The  total  contribution  of  any  given  vj  to  all  w’s  is  \iy,  i.e.,  2  =  P-j- 

/=o 
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M- 1_ 

•  The  sum  of  all  contributions  to  any  given  w,-  is  one;  i.e.,  ^  ry  =  1. 

im  o 


3.  Pyramid  Construction 

An  input  image  of  size  M  by  N  can  be  resampled  to  any  desired  output  size  K 
by  L  using  the  uniform  resampling  method  of  Section  1.1,  applying  first  a  resam¬ 
pling  of  the  rows  (which  are  N-vectors)  to  have  lengths  L,  and  then  resampling  the 
resulting  L  columns  (each  of  which  are  A/-vectors)  to  have  lengths  K.  The  result  is 
the  same  if  the  columns  are  first  resampled,  followed  by  a  resampling  along  the 
resulting  rows.  A  complete  pyramid  structure  can  be  built  by  specifying  the  desired 
sizes  of  all  levels  above  the  base  level. 

In  Figure  3,  we  present  an  example  of  a  three-level  one-dimensional  pyramid. 
Note  that  the  values  in  higher  levels  are  weighted  averages  of  bottom-level  pixels, 
and  that  the  supports  of  the  linear  weighting  functions  can  overlap  in  lower  levels. 
In  the  immediately  preceding  level  in  a  one-dimensional  pyramid,  with  resampling 
as  given  in  Section  1.1  or  1.2,  support  functions  for  adjacent  pixels  can  overlap  in  at 
most  one  pixel.  However,  the  supports  from  lower  levels  can  have  arbitrarily  large 
overlaps,  depending  on  the  structure  of  the  intermediate  levels.  Burt’s  overlapping 
pyramids  [2]  also  exhibit  overlapping  supports,  although  in  the  one-dimensional  ver¬ 
sion  of  his  scheme  (using  five-tap  filters) ,  adjacent  pixels  in  one  level  share  supports 
from  three  pixels  in  the  immediately  preceding  level  (see  Figure  4).  Successive  lev¬ 
els  have  larger  overlaps  in  the  base  level.  However,  resampling  using  either  the 
uniform  resampling  or  weighted  resampling  of  Section  1.1  or  1.2  yields  weights 
determined  by  the  sizes  of  the  levels,  and  are  not  adjustable  in  the  same  way  that 
the  taps  in  the  Burt  pyramid  can  be  modified. 


Figure  3.  A  comparison  of  a  4-3-2  one  dimensional  pyramid  with  a  4-2  pyramid 
structure.  The  leftmost  pyramid  shows  the  weights  in  a  4-3-2  structure,  while  the 
middle  pyramid  shows  the  effective  weights  from  the  bottom  level  to  the  top.  The 
rightmost  pyramid  gives  the  weights  for  a  uniform  4-2  resampling. 
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Figure  4.  The  Burt  5-tap  pyramid. 


Using  the  interpolation  formula  and  sampling  method  formulation  of  resam¬ 
pling,  however,  we  can  generalize  the  uniform  and  weighted  resampling  methods  to 
permit  larger  overlaps  in  the  support  functions  of  the  resampling  kernels.  Rather 
than  defining  <}><(•*)’ s  and  *|>f(jr)’s  as  in  Section  2.2,  we  instead  simply  assume  that 
the  interpolation  kernels  <j>i(*)  form  a  partition  of  unity  of  the  domain,  and  that  the 
sampling  kernels  <!><(*)  all  have  unit  mass.  By  this,  we  mean  that 

2  <►<(*)  -  1 
1-0 

and 

/*,(*)<&-  1,  for  i  =  0,  •••,  M. 

In  this  case,  resampling  is  still  done  by  the  formula  w,  =  2ryv;»  but  now  tbe  rij's 
are  given  by 

nj  =  ftyi(x)$j(x)dx. 

In  all  cases,  the  total  of  all  contributions  to  a  given  pixel  from  all  pixels  at  some 
given  lower  level  will  be  unity. 

The  pyramid  resampling  scheme  can  be  used  for  building  successively  smaller 
size  levels  (as  in  a  Gaussian  pyramid),  or  for  expanding  a  level  to  a  larger  image. 
By  combining  a  subsampling  operation  for  contraction  with  expansion  operations, 
the  resamplings  may  be  used  for  building  analogs  to  the  Laplacian  pyramid  data 
structure.  The  construction,  analogous  to  Burt’s  formulation,  proceeds  as  follows. 

We  specify  a  sequence  of  decreasing  sizes  ( N;  by  A/<)  for  levels 
/,  i  =  0,  1,  •  •  •  ,  €  of  a  pyramid  structure.  Let  Go  be  the  No  by  Afo  initial  image. 
We  construct  the  Laplacian  pyramid  as  a  sequence  of  images  L,  of  size  N,  by  Mt, 
i  =  0,  1,  •  •  t.  Recursively,  for  i  <  €,  Gj  +  i  is  obtained  from  G,  by  resampling 
from  an  N,  by  Mi  image  to  an  N,  +  i  by  Mi+\  image: 

Gj  + 1  =  ^(N|XMi)-(Ni+iXMl+1)  (Gj)  , 

where  7Z  is  the  resampling  operator.  We  assume,  for  the  time  being,  that  uniform 
resampling  will  be  used.  Then  L/  is  defined  as  the  difference  between  G,  and  an 
expansion  resampling  of  G,  +  r. 
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U  =  Gi  -  7Z(Nh.i  *  M/  +  l)  -  (J V,  x  M,)  (Gj  +  i)  . 

Finally,  at  the  top  level,  we  define  L(  =  G*. 

As  with  the  Burt  Laplacian  pyramid,  the  original  image  Go  can  be  recon¬ 
structed  exactly  from  the  L,-’s  using  the  (obvious)  formulas 

G(  =  L( 

Gi- 1  =  Li- 1  +  x  A/,)  X  W(_!)  (Gi)  . 


Each  level  L/  represents  an  approximate  band-pass  filtered  version  of  the  original 
image,  where  the  widths  of  the  bands  are  determined  by  the  choices  of  the  relative 
sizes  of  the  levels.  Figure  5  shows  such  a  custom-built  Laplacian  pyramid  structure. 

Suppose  that  we  wish  to  use  weighted  resampling  in  the  construction  of  the 
Laplacian  pyramid.  It  is  then  desirable  to  invert  (in  some  sense)  the  weighted 
resampling  that  was  used  in  the  contraction  from  level  G,-  to  G,  +  i  for  the  expansion 
operation  in  the  construction  of  L,-.  The  same  weighted  expansion  should  then  be 
used  in  the  reconstruction  of  G,-  using  the  levels  L,  +  i,  •  •  •  ,  L(.  A  weighted  expan¬ 
sion  resampling  should  be  used  (even  though  it  is  not  absolutely  required  for  exact 
reconstruction)  so  that  the  levels  L,-  maintain  their  approximate  bandpass  charac¬ 
teristics.  What  weights  should  be  used? 

Assume  that  the  data  (vo,  *  •  •  has  been  resampled  using  weights 

(no,  •  *  •  .P-at-i)  to  obtain  the  A/-vector  (vt>o,  •  •  •  Recall  that  the  resam¬ 

pling  is  linear: 

N-l_ 

W,  -  2  rijVj  . 

j- o 

The  matrix  coefficients  ry  represent  the  contribution  fractions  of  coefficient  vj  to 
coefficient  Wj.  We  now  define  the  “backward  weights”  (if,  i  =  0,  •  •  •  ,  M—\, 
according  to 


n-i  rtJ 

7  =  0  P-j 

M- 1 

It  is  easy  to  see  that  2)  =  so  that  the  (3,’s  can  serve  as  weights  for  resampling 

<=o 

the  Af-vector  (wq,  •  •  •  ,wm- i)  to  an  N-vector  (vo',  •  •  •  ,vn-\').  Each  p, 


represents  the  sum  of  proportions  of  v/s  that  contributed  to  w,-.  Accordingly,  the 
expansion  using  weights  (Po,  *  •  •  ,P m  -  i)  leads  to  a  vector  (vq',  ■  •  •  v^_i')  which 


best  approximates  (vq,  •  •  *  ,v^  -  i),  in  the  sense  that  there  is  no  horizontal  skew¬ 


ing. 


For  example,  uniform  resampling  of  an  N-vector  to  an  Af-vector  corresponds  to 
weighted  sampling  where  the  weights  are  all  equal  to  p  =  MIN.  The  backward 
weights  are  then  also  all  equal,  and  have  value  p-1  =  NIM.  Thus  uniform  resam¬ 
pling  for  contraction  is  paired  with  uniform  resampling  for  expansion. 


As  another  example,  we  showed  in  Figure  2  the  weighted  sampling  of  a  4- 
vector  to  a  3-vector  with  weights  (1,  1,  1/2,  1/2).  The  resulting  backward  weights 
are  (1,  1,  2),  so  that  the  weighted  backward  resampling  becomes 
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Figure  5a 


Custom-Made  Pyramid! 


Figure  5b 

Figure  5.  A  custom-built  two-dimensional  Laplacian  pyramid  with  level  sizes  as 
marked.  Levels  are  enlarged  to  full  size  for  visibility.  Figure  5a  show  the  Gaussian 
pyramid,  while  the  Figure  5b  contains  the  levels  of  the  Laplacian  pyramid.  The  ab¬ 
solute  value  of  the  data  is  displayed,  with  larger  values  shown  more  darkly. 
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as  depicted  by  Figures  6a  and  6b.  The  total  effect  of  the  contraction  followed  by  the 
expansion  is  that  the  higher  weighted  nodes  were  exactly  reconstructed,  while  the 
lower  weighted  nodes  were  blurred. 

4.  Adaptive  Resampling 

Weighted  resampling,  as  introduced  in  Section  1.2,  has  the  advantage  that 
regions  with  “interesting”  activity  can  be  resampled  more  finely  than  “uninterest¬ 
ing"  regions.  This  can  be  accomplished  (in  one  dimension)  by  giving  interesting 
pixels  larger  weights.  In  this  Section,  we  suggest  an  interest  operator  for  obtaining 
the  resampling  weights,  and  also  discuss  extensions  to  two  dimensions  and  irregular 
tessellation  grids. 

4.1.  One-dimensional  Adaptive  Pyramid 

We  suggest  an  interest  operator  based  on  the  local  “busyness"  of  the  data.  It 
has  been  observed  that  in  human  perception  a  line  with  higher  “busyness”  seems 
longer  than  a  straight  line  segment  [6j,  as  in  Figure  7.  Here,  we  will  use  a 
smoothed  absolute  value  of  the  Laplacian  of  the  data  to  measure  “busyness."  A 
similar  operator  has  been  suggested  for  representation  of  intensity  information  by 
retinal  receptive  fields  [7]. 
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Figure  6a.  Backward  weights  associated  with  the  weighted  resampling  of  Figure  2. 
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Figure  6b.  Calculation  of  the  resampling  weights  using  the  backward  weights. 


Specifically,  given  a  signal  vo.vi,  •  •  •  ,Vf/~  1,  we  compute  the  absolute  Lapla- 
cian  values: 

h  =  I  vi- !  -  2 v,  +  vi+1  |  ,  i  =  N-2, 

and  then  smooth  the  results  by  any  reasonable  local  averaging  operator,  extrapolate 
to  include  the  endpoints  t  =  0  and  i  =  N—  1,  and  normalize  to  obtain  the  values 
H<>  i  “  0,  •  ••  ,N-\.  One  way  to  accomplish  the  smoothing  is  to  define  an  itera¬ 
tive  weighted  smoothing  as  follows: 

*>o  j  =  lj,  lsjsjV-2, 

*>0,0  -  *1. 
bo,N-\  -  h-i  . 
and  recursively  set 


*>i+lj  **  1  +  \btJ  +  j*>,j  +  i  ,  1  S  S  N-2  , 

*>i+l,0  =  J*>», o  +  ^"*>f,  l  . 

*>»+l,N-l  =  J*>i,tf-2  +  J*>r,A>-l  • 

Finally,  set 

Nbr,i 

Hi  =  *17 

2  bTj 
J- o 

for  some  fixed  integer  T  >  0  representing  the  amount  of  smoothing.  The  m’s  ®re 
used  as  weights  for  resampling  the  v,’s. 

An  adaptive  custom-made  one-dimensional  Laplacian  pyramid  can  therefore  be 
easily  constructed.  The  sizes  of  the  levels  are  pre-specified.  The  signal  data  is 
given,  and  weights  are  obtained  from  the  busyness  measures  of  that  data. 
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Resampling  using  these  weights  results  in  the  data  immediately  above  the  base 
level.  Recomputing  new  busyness  measures  gives  new  weights  that  can  be  used  for 
the  next  level  resampling;  continuing  in  this  fashion  yields  a  Gaussian  pyramid  with 
adaptive  weights.  To  obtain  the  Laplacian  pyramid,  each  level  above  the  base  must 
be  expanded  to  the  size  of  the  previous,  level.  This  must  be  done  using  the  back¬ 
ward  weights,  as  described  in  Section  3.  The  resulting  expanded  levels  are  differ¬ 
enced  with  the  Gaussian  pyramid  data  at  that  level  to  obtain  the  Laplacian  data. 
Reconstruction  from  the  Laplacian  pyramid  structure  is  possible,  but  requires  that 
expansion  resampling  using  the  appropriate  backward  weights  at  each  level  be  used. 
This  means  that  the  complete  adaptive  Laplacian  pyramid  data  structure  consists  of 
the  Laplacian  data  at  each  level  together  with  the  backward  weights  needed  for  the 
expansion  resampling  at  each  level  above  the  base. 

4.2.  Two-dimensional  Adaptive  Pyramids 

The  two-dimensional  case  involves  technical  difficulties  not  present  in  the  one¬ 
dimensional  case.  The  problems  arise  due  to  the  fact  that  in  one  dimension,  weights 
can  be  converted  to  interval  lengths  to  decompose  the  domain,  whereas  in  two 
dimensions  it  is  not  clear  how  to  change  the  size  of  a  rectangular  region  in  relation 
to  a  weight  and  still  have  a  complete  decomposition  of  the  domain . 

One  solution  is  to  resample  in  each  dimension  separately.  However,  if  each 
row  or  column  is  processed  independently,  then  the  image  will  lose  its  structure  in 
the  sense  that  neither  connectivity  or  convexity  of  shapes  will  be  preserved.  Since 
the  weights  on  one  row  may  be  independent  and  different  from  the  weights  on  a 
neighboring  row,  nearby  pixels  may  give  principle  contributions  to  distant  pixels  in 
the  next  level.  Since  this  behavior  is  typically  unacceptable,  we  suggest  a  scheme 
where  all  rows  and  all  columns  use  the  same  weights.  In  this  case,  the  row  weights 
will  be  the  average  row  vector  of  a  busyness  matrix,  and  the  column  weights  will  be 
the  average  column  vector  of  the  same  matrix. 

Specifically,  we  compute  a  “busyness  measure”  fry  at  each  pixel  (ij)  in  the  N 
by  M  image.  Suppose  that  we  wish  to  resample  to  an  N'  by  M'  image.  Then  the 
row  weights,  used  to  resample  an  M-vt ctor  to  an  M '-vector,  are  given  by 

N- 1 

M'-'Zbij 
=  <-o 

N- 1  M- 1 

X  X  b* 

i- 0  k-0 

and  the  column  weights,  used  to  resample  each  column  N-vector  to  an  W '-vector, 
are  given  by 

Af-l 

N'-  X  bu 

j-0 

-  yrrferr  • 

X  X  bkj 

k-0  j-0 

The  same  weights  (p-o,  •  •  •  .p-m-i)  are  used  for  resampling  every  row,  and  the 
weights  (vq,  •  •  •  are  used  for  resampling  the  columns. 
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Once  again,  a  possible  method  for  computing  “busyness  values”  is  to  smooth 
absolute  values  of  the  Laplacian  of  the  image  data. 

Asa  trivial  example,  suppose  that  we  wish  to  resample  a  4  by  4  image  to  a  2 
by  2  image,  and  that  the  busyness  matrix  has  form 

0  0  0  0 
0  0  2  2 
0  0  2  2 
0  0  2  2 

Then  the  row  weights  are  given  by  (0,  0,  1,  1),  and  the  column  weights  are 
(0,  2/3,  2/3,  2/3).  The  row  resampling  will  then  work  as  shown  in  Figure  8a,  and 
the  column  resampling  is  as  shown  in  Figure  8b.  Note  that  the  result  is  that  the  3  by 
2  bottom  right  subimage  will  be  resampled  to  a  2  by  2  image,  and  that  the  other  pix¬ 
els  will  be  ignored. 

A  contrasting  example  arises  if  the  busyness  matrix  has  the  form 
2  2  0  0 
2  2  0  0 
0  0  2  2 
0  0  2  2 

In  this  case,  the  row  and  column  weights  are  both  (1/2,  1/2,  1/2,  1/2),  and  the 
resampling  for  both  the  rows  and  columns  will  be  uniform.  The  non-busy  quadrants 
(the  upper  left  and  lower  right  portions)  of  the  4  by  4  image  can  not  be  contracted 
without  distorting  shapes  in  the  remaining  portions  of  the  image. 

Figure  9  shows  an  adaptive  custom-made  pyramid  of  an  image.  Note  how  the 
interesting  regions  in  the  image  become  stretched  into  larger  windows  at  the  higher 
levels. 

5.  Irregular  Tessellations 

The  non-adaptive  versions  of  the  resampling  methods  given  in  previous  sec¬ 
tions  easily  extend  to  grid  arrays  other  than  rectangular  lattices.  Hexagonal  grids 


w0  W2 


Figure  8a.  Row  resampling  using  adaptive  weights  of  the  simple  example. 
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Figure  8b.  Column  resampling  using  adaptive  weights  of  the  simple  example 


Figure  9a 


Figure  9b 

Figure  9.  A  six  level  adaptive  pyramid,  with  level  sizes  as  marked.  Figure  9a 
shows  the  levels  of  the  Gaussian  adaptive  pyramid,  and  the  Figure  9b  shows  the  La- 
placian  adaptive  pyramid. 


are  often  suggested  as  a  useful  arrangement  for  image  processing.  There  are  hexag¬ 
onal  grids  where  the  cell  sizes  expand  with  increasing  distance  from  a  central  cell 
[8],  and  such  grids  may  have  a  basis  in  human  retinal  cell  distributions  [9].  We 
might  also  imagine  a  random  placement  of  sample  points,  with  local  density  of  sam¬ 
pling  points  prescribed  by  some  rule  (perhaps  adaptive  to  the  image  data).  We  will 
first  discuss  the  general  case,  and  then  focus  on  a  case  where  sample  point  locations 
are  described  on  a  polar  grid. 
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5.1.  General  case 

In  all  cases,  the  locations  of  the  sampling  points  are  associated  with  either  a 
regular  tessellation  or  irregular  tessellation  (a  decomposition)  of  the  image  domain 
into  disjoint  cells.  Each  sample  v*  represents  its  corresponding  cell  C and  the 
union  of  all  such  cells  covers  the  domain.  If  a  tessellation  is  not  provided  as  part  of 
the  sampling  locations,  then  one  can  be  provided  by  constructing  the  Voronoi 
decomposition  associated  with  the  points  [10]. 

Given  one  such  decomposition  of  the  domain,  say  {C  df-To1,  and  another 
(perhaps  coarser)  decomposition,  say  {PjJfLo1.  then  we  can  define  a  resampling  of 
data  defined  on  the  C  -grid  to  the  P-grid  as  follows. 

We  are  given  an  ^-vector  (v<j,  ■  •  •  of  data,  each  vj  representing  a  value 

on  cell  C  j.  We  define  resampled  data  (wo.  •  •  •  ,w^_j)  by 
N—  1 

w>  =  2  ruvj 
j~  i 

where 


rij  ■  f Styi(x,y)bj(x,y)dxdy  . 


Here  we  have  in  mind  interpolation  functions  »J»y(x,y)  that  are  characteristic  func¬ 
tions  of  the  corresponding  cells  C  j,  and  sampling  functions  that  are  characteristic 
functions  of  the  cells  P,  normalized  to  have  unit  mass.  That  is, 

<Ji;(*,y)  =  {  1  if  (**y)*Cj 
^0  otherwise  , 


and 


4>i(*.y) 


l/area(P<)  if  (x,y)«P< 
0  otherwise  . 


Area,  of  course,  is  measured  by  integration 


area  (P.)  =  ff  dxdy  . 


More  general  interpolation  and  sampling  kernels  can  be  envisioned,  but  the  simplest 
such  kernels,  the  characteristic  function  of  the  cells  as  described  here,  should  prove 
adequate  for  effective  resampling. 


5.2.  Polar  sampling 

We  now  discuss  a  special  sampling,  the  polar  sampling,  which  is  important  for 
snapshot  visual  perception.  In  polar  sampling,  the  sample  points  are  located  on  the 
intersections  between  a  set  of  rays  and  a  set  of  concentric  circles.  Some  models  of 
human  retinal  receptor  distribution  describe  samples  points  in  this  way,  where  reso¬ 
lution  falls  off  with  the  distance  (eccentricity)  from  the  fovea.  Sometimes  linear 
fall-off  of  the  sampling  rate  is  assumed,  which  leads  to  placement  of  grid  points  on 
concentric  circles  whose  radii  increase  as  log(A+M),  4*1,2.  ■  ,  where  A  and  h 

are  constants  [9].  Placement  of  the  circles  with  radii  increasing  exponentially  [11] 
and  also  simply  increasing  linearly  [12]  are  also  known  (see  Figure  10). 
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We  can  consider  resampling  a  polar  grid  by  reducing  the  number  of  concentric 
circles,  or  we  can  reduce  the  number  of  rays.  When  the  concentric  circles  are 
reduced,  we  simply  resample  the  grid  points  along  each  ray  independently  of  the 
other  rays.  Resampling  a  ray  is  straightforward  using  the  methods  of  Section  2. 
When  the  number  of  rays  is  reduced,  the  points  on  each  concentric  circle  can  be 
resampled  independently  of  the  other  circles.  To  do  this,  the  notion  of  one¬ 
dimensional  resampling  of  a  line  segment,  as  described  in  Section  2,  has  to  be 
extended  to  resampling  along  a  circle.  However,  the  methods  of  Section  2  extend 
easily,  using  distance  based  on  radian  measure  in  place  of  one-dimensional 
Euclidean  distance. 

If  both  the  number  of  concentric  circles  and  the  number  of  rays  are  to  be 
reduced,  we  can  then  do  the  resampling  first  by  reducing  the  number  of  circles,  and 
then  by  reducing  the  number  of  rays.  As  with  two-dimensional  resampling,  this 
process  is  unaffected  if  we  exchange  the  order  of  resampling.  Further,  the  optimal¬ 
ity  properties  of  the  weights  used  for  resampling,  alluded  to  in  Section  1,  extend  to 
polar  resampling.  Indeed,  because  the  polar  grid  coordinate  locations  can  be 
decomposed  into  a  cross  product  of  circles  and  rays,  a  separable  kind  of  adaptive 
resampling  is  possible  for  polar  grid  pyramids,  similar  to  the  adaptive  two- 
dimensional  rectangular  pyramids  discussed  in  Section  4.2. 

6.  Summary 

Using  the  anti-aliasing  method  common  in  computer  graphics  leads  to  an  idea 
for  one-dimensional  resampling  of  N  points  into  M  points.  Extending  this  idea  to 
two  dimensions  easily  establishes  a  method  for  building  pyramids  with  arbitrary 


Figure  10.  Polar  sampling,  where  sampling  points  are  on  the  intersection  points 
between  lines  and  circles.  The  left  sampling  is  with  exponential  radius  growth,  while 
on  the  right  the  growth  of  the  radius  is  linear. 
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sizes  specified  for  each  level.  We  thus  see  that  the  use  of  pyramids  with  dimensions 
given  by  powers  of  two  is  an  unnecessary  restriction  on  the  construction  of  the 
pyramids.  The  question  as  to  what  size  levels  are  most  appropriate  is  left 
unanswered,  and  depends  upon  the  application  and  empirical  experiences. 

We  have  also  investigated  the  idea  of  adaptive  resampling.  The  idea  is  easy  in 
one  dimension,  and  allows  for  arbitrary  nonlinear  stretching  along  the  line  to  be 
resampled.  We  gave  one  particular  busyness  measure  to  use  as  a  basis  for  deciding 
on  the  stretching.  In  two  dimensions,  things  are  more  complicated,  and  we 
compromised  by  permitting  only  a  “separable  stretching,”  where  the  rows  are 
resampled  using  one  set  of  weights,  yielding  the  same  stretching  for  all  rows,  and 
then  the  columns  are  resampled  using  a  single  set  of  weights. 

Finally,  we  note  that  the  idea  extends  to  irregular  tessellations,  where  sample 
points  can  be  randomly  placed,  or  placed  on  polar  grid  or  hexagonal  patterns. 
Ideally,  adaptive  resampling  would  allow  one  to  dynamically  place  additional  points 
in  regions  with  high  busyness,  but  our  method  does  not  easily  extend  to  the  case  of 
adaptive  placement  of  resample  point  locations  in  two  dimensions. 
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